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ON  THE  APPLICATION  OF  OPTICAL  SIGNAL  PROCESSING 


TO  THE  SOLUTION  OF  RICCATI  EQUATIONS 


I.  INTRODUCTION 

Riccati  equations  are  coupled,  nonlinear  differential  equations  which 
often  arise  in  the  analysis  and  design  of  modem  communication  and  control 
systems.  In  communications,  Riccati  equations  are  encountered  in  problems 
of  optimal  filtering  [l  J.  In  controls,  these  equations  arise  in  the  linear 
regulator  problem  [2],  in  servo-mechanism  design  and  Kalman  filtering  [3]. 

The  solutions  of  the  Riccati  equations  has  been  extensively  studied 
by  numerous  investigators.  The  general  solution  to  the  Riccati  equations 
is  quite  difficult,  because  these  equations  may  be  of  large  dimensions  (of 
order  50  or  more)  and  are  coupled,  nonlinear  equations.  The  advent  of 
large  scale  computers  has  permitted  some  progress  in  the  solution  of 
Riccati  equations,  but  the  computational  burden  is  still  quite  severe. 
Thus,  real-time  or  near  real-time  solution  of  these  equations  by  digital 
computations  is  presently  not  feasible  except  for  limited  cases  where  the 
dimension  of  the  equation  is  small.  Recently,  some  computer  algorithms 
have  been  presented  for  the  algebraic  Riccati  equation  (ARE),  which  is  the 
steady  state  solution.  Among  these  algorithms  are  those  proposed  by  Pappas 
and  Laub  [4J  using  the  generalized  eigenvalue  approach,  and  by  Laub  [5] 
where  he  uses  a  variant  of  the  classical  eigenvector  approach,  by  the 
so-called  Schur  vectors  method. 

Optical  signal  processing,  because  of  its  capacity  for  high  speed 
parallel  operations  using  relatively  simple  hardware,  has  emerged  as  a 
promising  candidate  for  solving  problems  in  modern  communications  and 
controls.  Various  forms  of  optical  signal  processing  architecture  have 
evolved  from  systems  designed  for  vector-matrix,  matrix-matrix  operations 
[6-9  j  and  eigenvector  determination  [lG,11J.  These  systems  are  often 
Manuscript  approved  June  13,  1983. 
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hybrid  electro-optical  systems  because  the  hardware  include  digital  and 
analog  electronics.  In  the  specific  area  of  modern  controls,  Casasent  et 
al  [l2J  described  an  iterative  optical  processor  (IOP)  for  solving  the  ARE 
by  implementing  a  revision  of  the  Kleinman  and  Richardson  algorithms.  This 
IOP  architecture  required  a  2-D  spatial  light  modulator  (SLM)  mask  with  a 
typical  refresh  time  of  30  msec,  thus  the  processing  speed  is  limited  by 
this  SLM  mask  refresh  time.  In  their  search  for  systems  with  higher 
speeds,  Casasent  et  al  described  an  architecture  [l3J  which  avoids  the  use 
of  2-E  SLM  masks.  This  system  was  applied  to  the  problem  of  the  discrete 
Kalman  filtering,  which  is  also  suitable  for  computer  solutions.  While  the 
discrete  Kalman  Filtering  is  useful  where  discrete  measurements  samples  are 
encountered,  continuous  Kalman  Filtering  is  important  when  a  continuous 
process  must  be  estimated  or  observed,  for  example,  in  the  shape  of  large 
space  structures. 

In  this  report,  we  describe  an  optical  signal  processing  technique 
for  the  solution  of  the  continuous  Riccati  equation  in  Kalman  filtering. 
Unlike  previous  algorithmic  techniques  [4,5,  and  1 2 J ,  the  optical  technique 
proposed  here  will  give  the  transient  as  well  as  the  steady-state 
solutions.  Using  the  great  power  of  parallel  processing  inherent  in 
optical  system,  the  solution  of  Riccati  equations  is  expected  to  be  in  near 
real-time. 


II.  RICCATI  EQUATION  IN  KALMAN  FILTERING 

In  Figure  1,  the  simplified  functional  block  diagram  of  a  continuous 
Kalman  filter  is  shown.  In  this  figure,  jc( t)  is  a  N-dimensional  vector 
which  represents  the  state  of  a  system  of  interest,  2(t)  is  the  estimate  of 
x(t),  which  is  computed  from  the  set  of  Nxl  measurements,  _z(t).  F(t)  and 
h(t)  are  NxN  matrices  which  are  called  the  system  and  measurement  matrices, 
respectively.  The  system  and  measurement  noise  vectors,  v(t)  and  v(t),  are 
Nxl  vectors.  In  the  right  side  of  figure  1,  K(t)  is  the  NxN  "Kalman  gain" 
matrix,  and  it  is  known  that  K(t)  is  dependent  upon  F(t)  by  [3J» 


2 


(1) 


K(t)  =  P(t)  HT(t)  R-1(t), 


T 

where  P(t)  is  the  solution  of  the  fiiccati  equation,  H  (t)  denotes  transpose 
of  H(t),  R_1(t)  is  the  inverse  of  B(t).  Id  (l)  we  assumed  v( t)  and  w(t) 
are  uncorrelated  and  R(t)  is  symetrical.  Thus  from  [3J 


P  *  FP  +  PFT  +  GQGT  -  PHT  R_1  HP,  P(tQ)  given,  (2) 

where  P  denotes  a  time  derivative,  P(t  )  is  the  known  initial  value  matrix 

o 

and  w(t)  and  v( t)  are  assumed  Gaussian  processes  with 
E[w(t)j  =  0,  E  [w2(t)]  =  Q(t) 

(3) 

E[v)t)j  =0,  E  [v2(t)]  =  R(t), 

with  E[  ]  denoting  taking  the  expected  value.  The  initial  condition  P(tQ) 
is  defined  by 


E  [(x(o)  -  x(o))(x(o)  -  x(o))T]  , 


1,4) 


where  E[x(o)j 


x(o).  The  correlation  between  w(t)  and  v(t)  is  given  by 


E  [w( t)  vT(t)]  =  C. 


(5) 


If  we  assume  R(t)  is  positive  definite,  then  R(t)  is  invertible  and  the 
solution  of  (1 )  is  assured. 
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III.  SOLUTION  OF  RICCATI  EQUATION 


The  system  of  N  nonlinear  coupled  equations  of  (2)  is  quite  difficult 
to  solve.  However,  if  the  coefficients  of  (2)  (i.e.,  F,  G,  H,  and  R)  are 
constants  (or  piecewise  constants),  then  by  considering  the  following 
transforms , 

X  =  Pjr  (6) 

and 

y  =  -FTy  ♦  kVkF  y  ,  (7) 

and  taking  a  derivative  of  (6)  with  respect  to  time  will  yield 

•  •  • 

X  =  py  +  py 

,  (pp  +  ppT  +  GCG1  -  PHTR_1HP)  y  +  P[-FTy  +  HTR“1HPy]  (8) 

=  FX  +  GQGTy. 

From  (7)  and  (6),  it  is  clear  the  matrix  equations  in  A_  and  £  can  be 
written  as 


which  is  a  2N  order  linear  system  of  equations,  that  is  much  simpler  than 

the  system  of  N  nonlinear  coupled  equations  of  (2).  From  linear  system 

theory,  it  is  known  the  system  of  (9)  can  be  written  as  *(t)  =  exp  (Et)  * 

$(t  +r,  t  ),  where  *(t)  is  the  transition  matrix.  Thus  we  can  write  [>]  as 
o  o 


5 


r  "i 

— 

—  — 

z(  V*) 

s 

V(,) 

Vx) 

*(t0} 

A(t  +t) 

0 

V"5 

*xx(,) 

X(t  ) 

0 

where  4(t)  is  partitioned  into  square  NxN  matrices.  Using  the  relation  in 
(6)  and  (10),  it  can  be  shown 


Mt0+-r)  =  p(tQ+-t)  jKt0+0* 
from  which 


(11) 


P(t0+-r) 


X(t  ) 

VT)  *  *>»(,)  ITT) 

xTTT 

•  (,)  +  a  (T)  ~yA- 

yy  yx  iTt07 


(12) 


P(Vt) 


+  P^o^  l-*yy(T)  +  V(T)  P(tO}j  " 


-1 


In  (12),  P(tQ)  is  assumed  known  a  priori  and  the  matrix  multiplication  and 
addition  may  be  performed  by  digital  techniques  or  by  combination  of 
digital-optical  techniques.  To  solve  equation  12  we  note  it  is  in  the  form 


P  *  VS 


(13) 


This  matrix  equation  can  be  solved  by  an  iterative  procedure  without  having 
to  invert  S.  Suppose  we  let  F  be  the  (i+l)th  iteration.  Then 

Pi+1  *  V  +  P.  (I-S),  (14) 
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where  I  is  the  identity  matrix  and  P.  is  the  i*'*1  iteration.  When  P.  .  =  P., 

1  i+I  l 

(14)  reduces  to  (13)  and  the  system’s  output  is  P.  Thus,  the  solution  P  is 
obtained  without  having  to  invert  S. 

IV.  DESCRIPTION  OF  OPTICAL  SYSTEM 


Consider  Fig.  2,  which  shows  the  proposed  optical  system  for  solving 

(13).  This  system  is  called  the  outer  product,  matrix-matrix  multiplier, 

iterative  optical  processor  (OIOP)  and  it  is  designed  to  solve  the  iterative 

algorithm  of  (14).  This  system  consists  of  the  outer  product  matrix-matrix 

multiplier  (due  to  Athale  and  Collins  [l4j)  which  computes  the  matrix 

product  P.(I-S),  an  electronic  adder  to  add  W  and  a  comparator  which 

compares  the  1  iteration  (P^)  with  the  (i+l)  iteration  (Pi+1). 

Functionally,  we  start  by  setting  P^*Pq  as  the  initial  value  of  P.  This 

value  P  is  then  summed  with  W  to  form 
o 


P 


i+1 


P  +  W. 
o 


(15) 


Assuming  W  f  0,  the  comparator  would  decide  P^+1  i  Pi  and  the  interation 

Pi+2  “  Pi+^(I-S)+W  would  be  performed.  This  iterative  procedure  continues 

until  P  =  P.  at  which  time  the  comparator  decides  the  solution  is 

1  C 

obtained.  We  now  proceed  to  discuss  functions  of  the  outer  product 

matrix-matrix  multiplier  and  the  comparator  circuit. 


A.  Outer  Product  Matrix-Matrix  Multiplier 


The  outer  product  matrix-matrix  multiplier  system  was  introduced  by 
Athale  and  Collins  [  1 4 j  to  perform  optical  matrix-matrix  multiplication. 
This  system  consists  of  two  I-B  acousto-optic  (AO)  cells  placed  orthogonal 
to  each  other  and  a  square  array  of  NxN  detectors,  where  N  is  dimension  of 
the  matrix.  This  system,  as  shown  in  Fig.  3,  operates  in  the  following 
manner.  Let  C*AE  be  the  matrix-matrix  multiplication.  For  the  case  where 


7 


DETECTOR 

ARRAY 


A  and  B  are  5x5  square  matrices,  the  ith  column  of  A  and  the  ith  row  of  B 
are  inserted  in  sequence  into  the  respective  AC  cells.  At  the  instant  when 
both  cells  are  filled,  the  light  source  S  is  strobed  to  effectively  freeze 
the  contents  of  the  AO  cells.  Thus  the  output  detector  will  realize  the 
outer  product 


b.,  b._  b .  b.,  b.= 
il  i2  i3  i4  i5 


1 1  il 

a. ,b.0 

1 1  i2 

a1 ibi3 

*tibU 

a  ,b.c 
1 1  x5 

„ .  b .  . 
2i  il 

a2i\2 

a2ibi3 

82ibi4 

a^.b.c 
£.1  l5 

3ibi1 

a3ibi2 

a_ ,b._ 
3i  i3 

a_.b.  . 
3i  i4 

a_ . b . c 
21  l5 

4ibi1 

a .  .b 

4i  i2 

a .  .b., 
4i  x3 

a  b.  . 

4i  i4 

a.  .b.c 
4i  i5 

5ibi1 

a  b.„ 
5i  i£ 

ac .  b._ 
5i  i3 

a_  .  b .  . 
51  i4 

ac . b . _ 
5i  i5 

If  successive  outer  products  between  the  columns  of  A  and  rows  of  B  (for  i 
=  1,2, 3, 4, 5)  are  formed  (as  shown  in  (16))  and  the  results  summed  at  the 
detector  array  by  using  time  integrating  detector  elements  (such  as  charge 
coupled  devices)  the  matrix-matrix  product  C  =  AB  will  be  realized.  To 
perform  rapid  matrix-matrix  multiplication,  the  row-column,  vector-vector 
multiplication  indicated  by  (16)  can  be  performed  by  inserting  the  matrix 
elements  in  parallel.  Thus  we  feed  the  columns  of  A  (a  )i*1 ,2, . .  .5) ,  ana 
the  rows  of  £  (b^i  *  1,2,...  5}  in  parallel  to  the  respective  AO  cells  (as 
shown  in  figure  4).  Vhen  a^'s  and  b^'s  are  properly  positioned  in  the  AC 
cells  the  light  source  is  strobed  to  provide  the  detector  array  output 
shown  in  (16).  Successive  column- row  outer  products  are  formed  and  summed 
at  the  detector  array,  as  previously  described,  to  obtain  the  matrix-matrix 
product  C  *  AB.  For  large  dimension  square  matrices  of  order  N  the 
parallel  feed  system  of  Fig.  4  is  N  times  faster  than  the  series  feed 
system  of  Fig.  3. 
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Eecause  the  OIOP  performs  operations  on  positive  real  information,  bipolar 
or  complex  data  must  be  first  transformed  or  rescaled  to  contain 
positive- real  data  components. 


B.  Bipolar  Lata 

Casasent  [l2j  shows  a  method  to  handle  bipolar  data.  Let  us  consider 
the  vector  matrix  product 


y  *  Hx 


(17) 


where  jr  and  x  are  Nxl  ve  era  and  H  is  a  NxK  matrix.  The  concept  can 
easily  be  extended  to  matrix-matrix  multiplications.  Let  h  be  an  element 
of  H.  We  define  a  new  matrix  B  with  positive  real  elements,  whereby 


b 

mn 


(16) 


where  h  and  h  are  respectively  the  minimum  and  maximum  values  of  H. 

clearly  CK  b  _<1 .  To  handle  the  bipolar  data,  we  decompose  each  element 
mn  ^  _ 

of  the  bipolar  input  vector  x  into  its  positive  x.  and  negative  x 
components  to  form 


a 


1  m 


C.5  (xm  ♦  |xj) 


(19a) 


a„  *  0,5  +  I ^ I ) t  (19b) 

cm  m  m 

where  m  *  1,2,...N.  From  ( 1 9 ) »  *  _x  and  are  both  positive  real 

and  it  can  be  noted  that 


x 


±1 


—2 


(20) 
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Equation  (19b)  above  differs  from  Eq.  6  of  Eef.  [l2j,  which  indicates  a  = 

2m 

0.5(x  —  | x  | )  that  is  apparently  in  error.  To  illustrate  the  use  of  (20) 

in  in  — 

by  a  simple  example,  consider^  «  (x1 ,  x2)  =  (5,-5).  Then  by  (19)  above 
an  *  C.5  (5  +  |5|)  =  5  ,  =  0.5  (-5  ♦  I -5 1 )  =  0 

(21) 

a21  -  0.5  (-5  +  |5|)  -  0  ,  a22  -  0.5  (5  ♦  15 1 )  -  5 

T  T  T  T 

Thus  a^  -  (5,0)  and  _a 2  =  (0,5)  are  both  positive  real  and 

T 

_x  =  -  a_2  *  (5,-5)  is  satisfied. 

An  implementation  for  the  system  ^  “  Hi  in  terms  of  _a^  and  _a 2  can  be  given. 
First,  note 

H  =  (h  -  h)  B  +  h,  (22) 

thus 

£  m  (h  -  h)  E  (_a^  -  a_2)  +  _h  x  (23) 

The  desired  system  for  bipolar  data  in  vector  matrix  multiplication  as 
indicated  by  (23)  is  shown  in  Fig.  5. 

C.  Complex  Data 

A  method  for  handling  complex  data  in  outer-product  matrix-matrix 
multiplicaton  was  shown  by  Taraseviich  et  al  [l5J«  In  this  method  the 
outer  product  of  complex  vectors  _a  and  b  with  positive  real  components  are 
performed  as  follows 


Fig.  5  -  Functional  block  liagra.  o£  bipolar  v.ctor-«ttl*  ™Ulpllar 
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In  the  above  Implementation  nine  times  as  many  detector  elements  are 
required  and  additional  electronic  circuitry  is  needed  tc  sum  the 
individual  outer  product  components  to  obtain  the  desired  result  of  (25c). 

Kore  economical  use  of  the  optical  system  space-bandwidth  product  can  be 

made  if  the  complex  data  is  represented  by  its  biased  real- imaginary  form 

T 

(1 6 J.  For  example,  the  outer  product  &  _t  is 
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Let  a  *  a  ♦  ia. 
—  ~r  i 


(27a) 


b  *  b 
—  — r 


(27b) 


T 

H  *  a  b  »  Mr  + 


(27c) 


+  j  (^i  — r 


♦ 


where  _a^  amd  are  respectively  the  real  and  complex  components  and  j*/-1 . 
To  perform  the  outer  product  in  an  acousto-optic  processor  requires 
positive  real  quantities.  Let  us  add  bias  vectors  \  and  Id  to  vectors  _a  and 
_b,  respectively.  Thus  the  outer  product  of  biased  vectors  becomes 


a  +  a 
— r  — 


a  +  a 

— i  — 


(£ +  iT: 


T  A  T  * 

(a  b  +  a  b  +  ab  +  a  b  ) 
— r— r  —  r~  — —  r  —  — 

Ca^B^  +  jiJb  +  ab*  +  ab  ) 


m  »m  ^  m  •*Am 

(a  b .  +  a  b1  +  a  b .  *  ab  ) 

— iP-i  — r—  —  —l  — 

m  A  m  A*ni 

(a^b£  +  a ^b1  +  _a  _b^  +  _ab  ) 


(26) 


To  obtain  the  desired  fo,  of  (27c),  appropriate  quantities  are  subtracted 
from  (28). 


C.  Accuracy 

Computational  accuracy  in  optical  processing  system  is  of  concern. 
Standard  acousto-optical  processing  procedures  is  limited  to  an  accuracy  of 
6-1 C  tits.  Recently,  techniques  for  improving  accuracies  in  optical 
computing  have  been  presented  by  Athale  et  al  [l7j.  In  this  technique, 
elements  of  matrices  are  represented  in  binary  form  and  binary 
multiplication  via  outer  product  is  combined  with  the  outer  product 
decomposition  of  matrix  multiplication  to  create  high  accuracy  optical 
matrix  multipliers.  To  multiply  two  KxK  matrices,  each  with  m-bit  dynamic 
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range  would  require  an  optical  system  with  an  mNxmN  space-bandwidth  product 
capability  and  a  dynamic  range  of  log2(mfJ)  bits  [  1 7 j - 

V.  SUMMARY  AND  DISCUSSIONS 

Vie  have  presented  an  outer  product,  matrix-matrix  multiplier  iterative 
optical  processor  (0I0P)  which  can  be  used  to  solve  the  matrix  equation 
P*WS  1  iteratively  without  having  to  invert  S  in  which  S  and  ¥  are  known 
matrices.  The  0I0P  can  be  used  to  solve  the  continuous  time  dependent 
Riccati  equation  that  arise  in  Kalman  filtering  problems  in  which  the 
coefficient  matrices  (i.e.,  F,  G,  H  and  R)  are  constants  or  piecewise 
constants.  A  functional  design  of'  the  0I0P  system  requires  a  matrix-matrix 
outer  product  multiplier,  electronic  summing  circuitry,  a  comparator 
circuitry,  and  a  feedback  loop.  The  iterative  algorithm  can  also,  of 
course,  be  implemented  digitally. 

Ve  have  also  briefly  discussed  methods  of  handling  bipolar  or  complex 
data  in  the  OIOP,  as  these  optical  processors  must  deal  with  non-negative 
real  data.  The  question  of  computational  accuracy  in  matrix-matrix 
multiplication  was  not  addressed  in  detail;  however,  the  technique  of 
binary  representation  using  outer  product  decomposition  [l7j  provide  a 
method  to  obtain  high  accuracy  at  the  expense  of  losing  some  of  processor 
space-bandwidth  product. 
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